install.packages(c("sf", "terra", "geodata", "rnaturalearth", "spdep",
"dplyr", "SpatialEpi", "wbstats", "flexdashboard",
"ggplot2", "viridis", "RColorBrewer", "patchwork", "DT",
"leaflet", "mapview", "leafpop", "leafsync", "rasterVis"))Geospatial Workshop PPSP 2024
Spatial Packages in R
Ensure the following R packages are installed beforehand.
install.packages("INLA",
repos = "https://inla.r-inla-download.org/R/stable", dep = TRUE)Load Packages
For this session, we will load some of the packages needed for plotting basic spatial data
library(sf)Warning: package 'sf' was built under R version 4.4.1
Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
library(mapview)Warning: package 'mapview' was built under R version 4.4.1
library(janitor)Warning: package 'janitor' was built under R version 4.4.1
Attaching package: 'janitor'
The following objects are masked from 'package:stats':
chisq.test, fisher.test
library(readxl)
library(dplyr)Warning: package 'dplyr' was built under R version 4.4.1
Attaching package: 'dplyr'
The following objects are masked from 'package:stats':
filter, lag
The following objects are masked from 'package:base':
intersect, setdiff, setequal, union
library(ggplot2)Warning: package 'ggplot2' was built under R version 4.4.1
library(DT)Warning: package 'DT' was built under R version 4.4.1
Making Maps with R
kelantan <- st_read("kelantan.shp")Reading layer `kelantan' from data source
`C:\Users\ACER\Desktop\geo_workshop\kelantan.shp' using driver `ESRI Shapefile'
Simple feature collection with 66 features and 6 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 371629.6 ymin: 503028.2 xmax: 519479.6 ymax: 690232.8
Projected CRS: Kertau (RSO) / RSO Malaya (m)
Plot 1
mapview(kelantan[,3])Plot 2
# Create a custom popup with multiple fields
popup_info <- paste0(
"<strong>MUKIM: </strong>", kelantan$MUKIM, "<br>",
"<strong>DAERAH: </strong>", kelantan$DAERAH, "<br>",
"<strong>Population: </strong>", kelantan$JUM_JANTIN
)
mapView(kelantan, zcol = "JUM_JANTIN", layer.name = "Population", popup = popup_info)Adding Spatial Attributes
Leptospirosis Data
lepto <- read_xlsx("leptospirosis.xlsx") %>% clean_names()
glimpse(lepto)Rows: 1,106
Columns: 21
$ diagnosis <chr> "Leptospirosis", "Leptospirosis", "Leptospirosi…
$ tahun_daftar <dbl> 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016,…
$ epid_daftar <dbl> 6, 9, 11, 12, 18, 18, 19, 17, 21, 23, 34, 34, 3…
$ age <dbl> 30, 23, 39, 43, 31, 34, 48, 55, 48, 39, 24, 45,…
$ alamat_semasa_kejadian <chr> "FELDA ARING", "LADANG U&I CIKU 8", "SYARIKAT S…
$ poskod <dbl> 18300, 18300, 18300, 18300, 18300, 18300, 18300…
$ latitude_rso <dbl> 478031, 459494, 441802, 488502, 496760, 441782,…
$ longitude_rso <dbl> 548141, 564966, 547551, 547701, 539072, 547566,…
$ latitude_wgs <dbl> 4.944824, 5.034372, 4.897434, 4.945540, 4.93648…
$ longitude_wgs <dbl> 102.3491, 102.1477, 101.9605, 102.3498, 102.454…
$ notifikasi_status <chr> "Daftar Kes", "Daftar Kes", "Daftar Kes", "Daft…
$ race <chr> "Foreigner", "Foreigner", "Foreigner", "Foreign…
$ kewarganegaraan <chr> "Bukan Warganegara", "Bukan Warganegara", "Buka…
$ gender <chr> "Male", "Male", "Male", "Male", "Male", "Male",…
$ nationality <chr> "INDONESIA", "INDONESIA", "INDONESIA", "INDONES…
$ klasifikasi_kes <chr> "Sporadic", "Sporadic", "Sporadic", "Sporadic",…
$ cara_pengesanan_kes <chr> "Pasif", "Pasif", "Pasif", "Pasif", "Pasif", "P…
$ jenis_rawatan <chr> "Wad Perubatan", "Jabatan Kecemasan & Kemalanga…
$ daerah <chr> "GUA MUSANG", "GUA MUSANG", "GUA MUSANG", "GUA …
$ mukim <chr> "CHIKU", "CHIKU", "GALAS", "CHIKU", "CHIKU", "B…
$ negeri <chr> "KELANTAN", "KELANTAN", "KELANTAN", "KELANTAN",…
Convert Leptospirosis data to spatial data
Use st_as_sf ( ) function to convert line listing data that contained Lat/Long attributes to spatial object
lepto <- st_as_sf(lepto,
coords = c("longitude_wgs", "latitude_wgs"),
crs = 4326)Confirm CRS is wgs84
st_crs(lepto)Coordinate Reference System:
User input: EPSG:4326
wkt:
GEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["EPSG",4326]]
Convert shapefile to RSO to match with Kelantan Map
lepto_2 <- st_transform(lepto, 3168)Plot map to see outlier
ggplot() +
geom_sf(data = lepto_2) +
ggtitle("Map of Leptospirosis") +
theme_bw()Plot the cases over the Kelantan map
overall_plot <- ggplot() +
geom_sf(data = kelantan) + #base map
geom_sf(data = lepto_2, color = "red", size = 0.5) + #cases in spatial data
ggtitle("Map of Leptospirosis Cases in Kelantan for 2016-2022") +
theme_bw() +
theme(plot.title = element_text(size = 12), strip.text = element_text(size = 12)) # cosmetic
overall_plot overall_plot +
facet_wrap(~tahun_daftar) + #to plot according to year
theme(axis.text.x = element_text(angle = 0, hjust = 0.5, vjust = 0.5)) +
theme(plot.title = element_text(size = 12), strip.text = element_text(size = 12)) +
ggtitle("Map of Leptospirosis Cases in Kelantan for 2016-2022")Using OSM as basemap
library(leaflet)Warning: package 'leaflet' was built under R version 4.4.1
leaflet(lepto) |>
addTiles() |> # Add default OpenStreetMap tiles
addCircleMarkers(
radius = 2, # Adjust marker size
color = "red", # Marker color
stroke = FALSE, fillOpacity = 0.8)Aggregated Data
For this exercise, we focus our analysis on Leptospirosis cases reported in 2016 only.
lepto_16 <- lepto_2 %>% filter(tahun_daftar == "2016")Joint point to polygon
st_crs(lepto_16)Coordinate Reference System:
User input: EPSG:3168
wkt:
PROJCRS["Kertau (RSO) / RSO Malaya (m)",
BASEGEOGCRS["Kertau (RSO)",
DATUM["Kertau (RSO)",
ELLIPSOID["Everest 1830 (RSO 1969)",6377295.664,300.8017,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4751]],
CONVERSION["Rectified Skew Orthomorphic Malaya Grid (metre)",
METHOD["Hotine Oblique Mercator (variant A)",
ID["EPSG",9812]],
PARAMETER["Latitude of projection centre",4,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8811]],
PARAMETER["Longitude of projection centre",102.25,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8812]],
PARAMETER["Azimuth of initial line",323.0257905,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8813]],
PARAMETER["Angle from Rectified to Skew Grid",323.130102361111,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8814]],
PARAMETER["Scale factor on initial line",0.99984,
SCALEUNIT["unity",1],
ID["EPSG",8815]],
PARAMETER["False easting",804670.24,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Engineering survey, topographic mapping."],
AREA["Malaysia - West Malaysia onshore."],
BBOX[1.21,99.59,6.72,104.6]],
ID["EPSG",3168]]
Count all leptospirosis cases for each mukim in 2016
count_lepto_sub_yr <- lepto_16 %>%
count(daerah, mukim) %>%
ungroup()Joining the count data to base-map of Kelantan
kelantan_count_lepto <- st_join(kelantan, count_lepto_sub_yr)Map the count of cases according to subdistrict in Kelantan
# Create a custom popup with multiple fields
popup_info <- paste0(
"<strong>MUKIM: </strong>", kelantan_count_lepto$MUKIM, "<br>",
"<strong>DAERAH: </strong>", kelantan_count_lepto$DAERAH, "<br>",
"<strong>Leptospirosis Cases: </strong>", kelantan_count_lepto$n
)
# Display the map with the custom popup
mapView(kelantan_count_lepto, zcol = "n", layer.name = "Count of Leptospirosis cases 2016", popup = popup_info)